sir <- function(t, y, pars){
# State variables
S <- y[1]
I <- y[2]
R <- y[3]
# Parameter values
beta <- pars["beta"]
gamma <- pars["gamma"]
# Equations
dS <- - beta*S*I
dI <- beta*S*I - gamma*I
dR <- gamma*I
out <- c(dS, dI, dR)
# Return list of gradients
list(out)
}
times <- seq(0, 0.25, by = 1/365)
paras <- c(beta = 200, gamma = 365/7)
init <- c(S = 0.999, I = 0.001, R = 0)
data_out <- ode(y = init, times = times, func = sir, parms = paras)
data_out <- as_tibble(data_out)
plot <- data_out %>% pivot_longer(cols = S:R,
names_to = "Variable", values_to = "Proportion") %>%
mutate(timeD =time*365) %>%
ggplot(aes(x = timeD, y = Proportion, col = Variable)) +
geom_line(size = 1) + theme_bw() + xlab("Time (days)")
plot
seir <- function(t, y, pars){
# Variables
S <- y[1]
E <- y[2]
I <- y[3]
R <- y[4]
#parameter
beta <- pars["beta"]
gamma <- pars["gamma"]
sigma <- pars["sigma"]
#equations
dS <- - beta*S*I
dE <- beta*S*I - sigma*E
dI <- sigma*E - gamma*I
dR <- gamma*I
out1 <- c(dS, dE, dI, dR)
list(out1)
}
#parameters
times1<-seq(0, 0.25, by = 1/365)
paras1 <- c(beta = 200, gamma = 365/7, sigma = 365/5)
init1 <- c(S = 0.999, E = 0, I = 0.001, R = 0)
#run code
data_out1 <- ode(y = init1, times = times1, func = seir, parms = paras1)
data_out1 <- as_tibble(data_out1)
head(round(data_out1, digits = 2))
## # A tibble: 6 × 5
## time S E I R
## <deSolve> <deSolve> <deSolve> <deSolve> <deSolve>
## 1 0.00 1 0 0 0
## 2 0.00 1 0 0 0
## 3 0.01 1 0 0 0
## 4 0.01 1 0 0 0
## 5 0.01 1 0 0 0
## 6 0.01 1 0 0 0
tail(round(data_out1, digits = 2))
## # A tibble: 6 × 5
## time S E I R
## <deSolve> <deSolve> <deSolve> <deSolve> <deSolve>
## 1 0.24 0.02 0 0.01 0.96
## 2 0.24 0.02 0 0.01 0.96
## 3 0.24 0.02 0 0.01 0.96
## 4 0.24 0.02 0 0.01 0.97
## 5 0.25 0.02 0 0.01 0.97
## 6 0.25 0.02 0 0.01 0.97
#plot
plot1 <- data_out1 %>% pivot_longer(col = S:R, names_to="Variable", values_to = "Proportion") %>%
mutate(timeD = time*365) %>%
ggplot(aes(x = timeD, y = Proportion, col= Variable)) + geom_line(size=1) + theme_bw() + xlab("Time(days)")
plot1
In comparison to the SIR model, the SEIR model (a) the peak of the infection is reached later than the SIR model. The peak is reached before the 25th day in the SIR versus the SEIR is is reached after the 25th day possibly even doubled. In comparison to proportion of susceptible once the out break goes extincted, they are both similar in that the proportion of susceptible will never reach 0 but will remain slightly above.
si <- function(t, y, pars){
# Variables
S <- y[1]
I <- y[2]
#parameter
beta <- pars["beta"]
gamma <- pars["gamma"]
#equations
dS <- - beta*S*I
dI <- beta*S*I
out2 <- c(dS, dI)
list(out2)
}
#parameters
times2<-seq(0,0.25, by = 1/365)
paras2 <- c(beta = 200, gamma = 365/7)
init2 <- c(S = 0.999, I = 0.001)
#run code
data_out2 <- ode(y = init2, times = times2, func = si, parms = paras2)
data_out2 <- as_tibble(data_out2)
head(round(data_out2, digits = 2))
## # A tibble: 6 × 3
## time S I
## <deSolve> <deSolve> <deSolve>
## 1 0.00 1.00 0.00
## 2 0.00 1.00 0.00
## 3 0.01 1.00 0.00
## 4 0.01 0.99 0.01
## 5 0.01 0.99 0.01
## 6 0.01 0.98 0.02
tail(round(data_out2, digits = 2))
## # A tibble: 6 × 3
## time S I
## <deSolve> <deSolve> <deSolve>
## 1 0.24 0 1
## 2 0.24 0 1
## 3 0.24 0 1
## 4 0.24 0 1
## 5 0.25 0 1
## 6 0.25 0 1
#plot
plot2 <- data_out2 %>% pivot_longer(col = S:I, names_to="Variable", values_to = "Proportion") %>%
mutate(timeD = time*365) %>%
ggplot(aes(x = timeD, y = Proportion, col= Variable)) + geom_line(size=1) + theme_bw() + xlab("Time(days)")
plot2
As we can see in the SI model, we are only using two parameters susceptible and infected individuals. Looking at the plot it shows that peak infection, intersection of the plot, is at around 20-21 days. because S and I are inverses of each other, as one loses the other gains, the differential equations then now equal 0.
sirs <- function(t, y, pars){
# State variables
S <- y[1]
I <- y[2]
R <- y[3]
# Parameter values
beta <- pars["beta"]
gamma <- pars["gamma"]
omega <- pars["omega"]
# Equations
dS <- - beta*S*I + omega*R
dI <- beta*S*I - gamma*I
dR <- gamma*I - omega*R
out <- c(dS, dI, dR)
# Return list of gradients
list(out)
}
times3 <- seq(0, 5, by = 1/365)
paras3 <- c(beta = 200, gamma = 365/7, omega = 365/30 )
init3 <- c(S = 0.999, I = 0.001, R = 0)
data_out3 <- ode(y = init3, times = times3, func = sirs, parms = paras3)
data_out3 <- as_tibble(data_out3)
plot3 <- data_out3 %>% pivot_longer(cols = S:R,
names_to = "Variable", values_to = "Proportion") %>%
mutate(timeD =time*365) %>%
ggplot(aes(x = timeD, y = Proportion, col = Variable)) +
geom_line(size = 1) + theme_bw() + xlab("Time (days)")
plot3
plot4 <- data_out3 %>% pivot_longer(cols = S:R,
names_to = "Variable", values_to = "Proportion") %>%
mutate(timeD =time*365) %>%
ggplot(aes(x = timeD, y = Proportion, col = Variable)) +
geom_line(size = 1) + theme_bw() + xlab("Time (days)") + xlim(0,125)
plot4
I ran my plot at the 5 year mark. I have added to plots to see it at the 5 year span and when we start to see the consistency between S,I and R. as can be seen in the plot. Because this in SIRS model, we can see how it already differs from the SIR model. In the SIR model when can see that R will continue to get closer to 1 while S and I get closer to 0. Where as in the SIRS model S,I, and R become constant at 0.25, 0.13, and 0.62 respectively, this is because individuals can become susceptible the infection will stay within the population.
Extra Credit: